%{
                                                     nearestNeighbors.m V1
                                                     Written by Rajeev Kant
                                                        Sept. 4th, 2020
                                    Coulombe Lab for Cardiovascular Regenerative Engineering
                            Center for Biomedical Engineering, Brown Univeristy, Providence, RI, 02903

This function calculates the nearest neighbors for a certain particle
within a defined search radius using a basic for-loop methodology. Note
that denser images with larger search parameters exponentially increase
run-time.

Inputs:
1. stats: a regionprops stats input with Centroid and PixelList fields from a binarized image of nuclei
2. path: file path
3. im: binarized nuclei image (in this case a subROI) corresponding to
regionprops stats
4. radius: search radius in %%%%%microns%%%%% to find neighboring particles
5. scale: Scale facotr
6. nameOut: naming convention for saved image

Outputs:
1. stats: Saves nearestNeighbor outputs into the input stats file
consisting of:
 1a. nnCount: Number of neighbors within a specified radius
 1b. nnCoords: Paired coordinates from particle of interest to nearest
 neighbors for graphing
 1c. nearestDist_um: Nearest distance in um to next closest particle. NaN
 for particles greater than the set radius.
2. Plot of segmented nuclei with lines connecting them to nearest neighbors


Original publication for citation: 
Kant, R. J., Bare, C. F., and Coulombe, K. L. K. "Tissues with 
patterned vessels or protein release induce vascular chemotaxis
in an in vitro platform", Tissue Eng. Part A, 2021.
doi:10.1089/ten.TEA.2020.0269

%}
function [stats] = nearestNeighbors(stats, path, im, radius, scale, nameOut)
%% Find nearest neighbors within specified radius

rad_px = radius*scale; % Converts input radius back into pixels for analysis

% Determine max array size for preallocation to minimize array size and optimize runtime
maxSize = 0;
for i=1:length(stats)
    len = length(stats(i).PixelList);
    if len>maxSize
        maxSize=len;
    end
end

% Initialize Variables
nnDists = zeros(maxSize, 1, length(stats)); % Distances of nearest neighbors
nnCoords = zeros(maxSize, 4, length(stats)); % Corresponding coordinates of nearest neighbors
nnCount = zeros(1, 1, length(stats)); % Total number of nearest neighbors within specified radius

% Compare each pair of particles and for those within specified radius,
% Compare individual pixels to find shortest distance between them
for i = 1:length(stats) % Iterate through each particle
    % Reinitialize placeholder values each loop
    partMins = ones(maxSize,1)*rad_px; % Set default values to rad
    minC = zeros(maxSize,4);
    nn = 0;
    for j = 1:length(stats) % Iterate through jth particle to be compared to ith particle
        if isequal(i,j) % Skip comparison if comparing particle to self
            
        else
            % Check centroid distance and remove particles that are 2*rad away to reduce analysis time
            cen = cat(1, stats(i).Centroid, stats(j).Centroid);
            cenDist = pdist(cen);
            if cenDist > rad_px * 2
                
            else % Calculate distance between each pair of points in the pixel lists of particle i and j and saves the absolute minimum distance and corresponding points
                for k = 1:size(stats(i).PixelList,1) % Iterate through each pixel in i's pixel list
                    for m = 1:size(stats(j).PixelList,1) % Iterate through each pixel in j's pixel list
                        X = cat(1, stats(i).PixelList(k,:), stats(j).PixelList(m,:));
                        dist = pdist(X); % Calculate distance between paired points
                        if dist < rad_px && dist < partMins(nn+1,1)
                            % Save and overwrite if next pair of points is closer than previous (default is rad) - nn+1 to reference index
                            partMins(nn+1,1) = dist;
                            minC(nn+1,1:2) = stats(i).PixelList(k,:); % Saves coordinates
                            minC(nn+1,3:4) = stats(j).PixelList(m,:);
                        end
                    end
                end
                if partMins(nn+1,1) < rad_px % Only increment nn if a value less than 10 was read in - filters out values rad > dist < 2*rad
                    nn = nn+1; % Increments number of neighbors counter after comparing pixel lists(k and m) beween ith and jth particle
                end
            end
        end
    end
    % Write minimums to 3D arrays prior to iterating to next particle
    partMins(nn+1:end,1) = 0; % Changes all other non NN counts from 10 to 0 (empty cells)
    nnDists(:,:,i) = partMins(:,:);
    nnCoords(:,:,i) = minC(:,:);
    nnCount(:,:,i) = nn;
end

%% Draw plot over im with lines connecting neighbors

iptsetpref('ImshowBorder','tight');
figure; imshow(im);
hold on;
for i=1:size(nnCoords,3)
    for j=1:nnCount(:,:,i)
        xc = [nnCoords(j,1, i) nnCoords(j,3, i)];
        yc = [nnCoords(j,2, i) nnCoords(j,4, i)];
        pl = plot(xc, yc);
    end
end

% Save image
saveas(gcf,fullfile(path, strcat(nameOut,'.tif')),'tiff');
close all;

%% Update stats and prepare summary output

% Read nnCount into stats
cellCount = num2cell(nnCount);
[stats(:).nnCount] = cellCount{:};

% Cut nnCoords and nnDists at nnCount (to remove excess 0's), calculate nearestDist and read all into stats
for i = 1:length(stats)
    n = nnCount(:,:,i);
    nC = num2cell(nnCoords(1:n,:,i));
    nD = nnDists(1:n,:,i)/scale; % Output for nnDists and nearestDist is converted back to �m
    nDCell = num2cell(nD);
    
    stats(i).nnCoords = nC;
    stats(i).nnDists_um = nDCell;
        
    if n == 0
        stats(i).nearestDist_um = NaN; % Read in nearestDist for particles with no neighbors as NaN
    else
        nDmin = min(nD);
        stats(i).nearestDist_um = nDmin;
    end
end

% End of nearestNeighbors.m